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ABSTRACT 
This thesis proposes a method of identifying the dynamics of ship 
angular motion at sea as the basis for ship motion estimation. Various 
sources of information are discussed. A digital model to simulate 
ship's motion is developed. Identification algorithms are presented 
with the results of their stochastic digital simulation. Sensitivity of 
identification error to sample rate was investigated. A plan for ship- 


board implementation utilizing a digital computer is presented. 
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INTRODUCTION 

The sophistication of modern shipboard systems is imposing 
accuracy requirements in meaSuring a Ship's attitude, A 10 mil error 
in yaw correction to a gun platform is a 200 yard error at 10 nautical 
miles. System time delays and additive measurement noiSe are key 
sources of system error, Automatic carrier landing systems require 
that the predicted touch down point on the deck be known, Present 
system effectiveness degenerates rapidly as sea states increase, 
When the wave-off decision point is reached by the aircraft, the auto- 
mated landing system must know the aircraft position relative to the 
desired impact point at the time of touchdown. If pitch and yaw can 
accurately be predicted, a true all weather capability may be antici- 
pated. For these and many other applications the accurate knowledge 
of present position and the prediction of future positions are key 
elements of the system environment. 

An accurate prediction of ship motion must necessarily involve 
information from ‘the forcing function. Proper interpretation should lead 
to an identification of the current dynamic structure of the ocean waves. 
This could yield a more accurate input to the prediction of sea states 
and ocean weather propagation. 

The French have done considerable work in the field of spectro- 
angular ocean wave analysis [1, 2 ]. However their work is success- 
ful only in predicting mean and extreme wave amplitudes and periods. 
The primary available inputs are local surface wind conditions dis- 
tributed over the ocean area. Workin this area is not yet sufficiently 
advanced to be applied to a real time estimation problem, 

The investigation of the statistical properties of ocean waves has 
been a well traveled route. R. L. Wiegel | 3 ] from the University 
of California has done extensive work using both surface and bottom 
moored pressure senSing units. Data reduction again shows normal 


distributions of wave period and amplitude in open waters. 


Ship motion estimation may be approached as a stochastic con- 
trol problem. Thus the ship is the plant and the ocean wave is the 
forcing function with normal Gaussian distribution. Pitch, roll, 
and yaw will be treated as separate, uncoupled subsystems dis- 
cretely sampled to facilitate digital simulation and data processing. 
Only the roll mode will be referred to in the remainder of this paper 
on the assumption that the other modes may be estimated by the same 


techniques. 


PLANT SIMULATION 
To evaluate a technique for ship motion estimation, it is first 
necessary to have a source of data. A fourth order plant composed of 
two cascaded identical second order systems was selected for motion 
simulation. This plant has known natural frequency ( w ) and damping 


coefficient (C ), 


] l 
G(s) = ia _ (1) 
s +20CW st+w s +2Cw st+w 
n n = hh n 

The plant may be represented by 

X = AX+BW_ where (2) 

0 l 0 0 0 
A. 0 0 1 0 ae = 0 

0 0 0 il 0 

2 
—W -4 a Peace 4.5) -4 lw l 
Lan n n n| 


Selection of © and Ww determines the dominant characteristics. The 
discrete recursive state equation in matrix form is 


X(k + 1) = ¢ X(k) + IF Wk). (3) 


W(k) is the discrete, white Gaussian excitation with zero mean. The 
desired system output is the discrete time sequence of roll angle 
corresponding to x) where 
X, (k) = HX(k), H=(100...). 
Additive measurement noise, v(k), may be included at the output to 
complete the simulation model, fig. 1. V(k) is assumed to be white 
Gaussian noise, uncorrelated with W(k). The process is described by 
z(k) = HX(k) + v(k). (4) 
It will be convenient to transform ¢ginto the canonic companion 


matrix form, ¢ *, using Lee's observability criteria [4 ], giving 


H ¢ 
D= y, 
Hf glee 
t 
1 0} I 
$*=Dg De = |--t--------== ) (5) 
ay an - a 


and the companion distribution matrix 
|e — a (6) 


From the work of Lee and Ho [5 J} this defines a new recursion equation 


with the same output values as equation (4) as seen below 


Y(k+1) = ¢@ * Y(k) + TI * W(k) 
zk) =H Y(k) + v(k) (7) 


where W(k), H, v(k), and z(k) are the same but the Y elements are 
different. 

Representing the simulation model in this form reduces both the 
number of computations and the time required to generate a sample set 
of data ona digital computer. Since the data set is unchanged by the 
foregoing manipulations, the dominant characteristics are preserved. 
Values of C= 0.7 and Ww = 0.2 7 were used as reasonable estimates 
for the roll dynamics of a moderate size ship. The excitation, W/(k), 
was a random number generating subroutine with zero mean and stan- 
dard deviation of one. Multiplying the random numbers by a con- 
stant permitted simulation of any desired sea state. An example of 
roll angles generated using a constant of twenty and a sample period 


of 0.1 second is shown in fig. 2. 
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Fig. 1. Block Diagram of the Simulation Model. 
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Fig. 2. Roll Angles vs Time (10 sec/in) Generated by the Simulation 
Model Excited with Discrete Gaussian Noise, o= 20, at the 
Sampling Interval, 0.1 sec. Dominant Period = 10 sec, 
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IDENTIFICATION 

The problem of identification arises when the set of differential 
equations describing the system dynamics are not known, uncertain, 
or are subject to periodic or continuous change. Ona large time 
scale the ship dynamics may be expected to change as fuel, Stores, 
ammunition, etc. vary in quantity and location. Short time changes 
may be anticipated due to changes in ships heading and speed. 

In state vector notation the recursion equation for ship roll 


motion (uncoupled) will have the form; 
X(k+1) = @X(k) + IT W(k) (8) 


where W(k) is the discrete Gaussian excitation. © is the fixed 

state transition matrix of rank N, the order of the system. The term 
fixed is used here on the assumption that the changes in dynamics 
mentioned above occur at such a sufficiently slow rate that @ may be 
adequately represented as time invariant over a finite time interval. 
Thus © must be periodically identified to follow ship motion within 


accuracy specifications. 
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IDENTIFICATION OF THE FREE DYNAMIC SYSTEM 


A free dynamic system may be represented by the equation 
Y(k+1) = @ * Y(k) . (9) 


Lee developed an identification scheme for a plant having no numerator 


dynamics [4 ]. Lee's development gave the following equations; 


oe a ae (10) 
where 
ee, oem “2 73 °° area) 
ames “3 «74 
SSL : ‘ aA : 
an ntl *On-1 eae ae On 


and z(k) = HY(k), a scalar. 
Using the same data, the order of the system may also be iden- 
tified. Choose some M greater than N and build an M xX L matrix, A. 


i I 2 Cie hs Sk eae 


i 2 L 
“9 Oe a a es 
\— 
a 9 9 lial tag pa) 


The product ata will be positive definite for L less than or equal to N 
and singular for L greater than N. Hence the order of the system is 
the maximum nonsingular value of L. 

A problem arises from the possible occurrence of numerator dy- 
namics which leads to residues (errors) in the identification. In z- 
transform theory the nth row of ¢ * corresponds to the denominator 


coefficients of the system z-transform. The numerator remains 
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unidentified. Hence, if there is an unknown non-scalar numerator in 
the real plant corresponding to a nonzero initial condition, errors will 
result. 

Identification using this method was investigated with excellent 
results by Setting an initial condition for roll angle and then releasing 
the plant. If excitation is applied while measurements are being taken, 
the plant is no longer a free dynamic system, Identification attempted 


with random excitation applied was inaccurate and unreliable. 
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IDENTIFICATION FOR THE STOCHASTIC CASE 
Using the transformed recursion equation in the new state space, 
Lt. Ralph Hudson's unpublished doctoral notes show the development 


of a straight forward identification scheme as follows: 


Y¥(k+1) = @ * Y(k) + wk) (11) 
where 
* ntl 
Y(k) = . Pill) ele), 
k-1 


and the scalar, W(k), represents the white Gaussian excitation. If 


equation (11) is post multiplied by Y(k)': 
¥(k+1)¥(k)! = @ *¥(k)Y(k)? + w (k)¥(k)! (12) 
Taking the expectation of both sides; 
EL Y(k+1)¥(k)"] = ¢ *EL Y(k)Y(k) J +ELw(k)¥(k) ‘1 (13) 


Notice that the EL Y(k+1)¥(k)* J is just the autocorrelation function, 
R( +), for tau equal to one. Similarly, EL Y(k)¥(k)*] is the auto- 
correlation for tau equal to zero. Therefore 3 x= R(1)R(0) - 

If the order of the identification is M and M <N,, the true system 
order, the inverse will exist. This is implied by the fact that an nth 
order linear system may be defined by N linear independent differ- 
ential equations. Hence any M <N of these equations are also in- 
dependent. If M >N independence is lost and R(0) of rank M >N is 
Singular. Thus if N is unknown, it may be identified by testing for 
the largest non-singular rank of R(0). 

For application to discrete linear (Kalman) filtering it is also 
useful to note that the © required for determining the gain matrix may 


be concurrently identified as follows; 


O=TF* EL wk)wk) Ir *) =EL wk) wk)‘ (14) 


Es 


Hence 


Q=E[ {v(k+l) - o*¥(k)} {¥(k+1) - 9 *¥(K)} 77 


which simplifies to 


Q=R(0) - ¢ *R(1). (15) 


This scheme was programmed in Fortran 63 for simulation on the 
CDC 1604 digital computer using data from the simulation model, 
Batch processing and recursive identification were investigated. 

1, Batch processing: A sample set of roll angles was generated 
and stored. The stored data was then batch processed to form the auto- 
correlation functions, R(0) and R(1), for the set. A matrix inversion 
routine was used to invert R(0). ¢@* was then identified as the matrix 
product, R(1)R(0) + 

If I’ * was restricted to a single non zero element, identifi- 
cation of the fourth order plant was good using 50 samples and ex- 
cellent for 500 samples. However, if I * had all non zero elements, 
residues due to numerator dynamics caused errors as great as 110% 
of the true values of the elements. Increasing the sample size from 
500 to 900 showed no improvement. 

2. Recursive identification: Hudson's identification scheme may 
be manipulated into the following recursion equations by uSing a matrix 


inversion lemma. 


o*(k+1) = ¢*(k) + (ve) >oa¥0e + 1) oa (ves) =% *()Y(K)) Y(k)  P (k) 
P(k+1) = P(k) (: eer vie Tl vie) vee 17>) 


where P(k) = Ro 


Let B (k) = (via Pam ¥(K) + 1) ay a scalar. 


Then 
b*(k+1) = $*(k) +B (k) (vee ‘ 3+(K)¥(k) | ¥(k)* P(k) (16) 


IG 


and 


P(k+1) = P(k) (1 - B (K)¥(4)¥(4)"P(K) ) ; (a7) 


The above relationships identify an N by N ¢ * matrix of Nn? elements. 
For the scalar observable case N(N-1) of these elements are already 
known to be zeros or ones. Hence N’ elements are identified to learn 
N unknowns in the nih row. Thus a more efficient method of the fore- 
going has been derived by Yu-Chi Ho | 5 ] and R.C.K. Lee [ 4,5 ] 


which recursively estimates the elements of the Nee row. 


(k+l) = (k) + P(k)¥(k) (204) 2 $ tk) "v¢k) ) B (k) (18) 
P(k) = P(k+1) (1 - B (k-1)¥(k-1)¥(k-1)Pk-1)) (19) 


where z(k+1) = the scalar measurement at stage (k+1) 


A A 
@(k) = (Nx1l) column such that ¢ (k)? is equal to the nt 


h 
row of ¢ *(k) 

Initialization of the recursive equations might be accomplished for 

A 
* = 

6 *(0) anes 1 F 
measurements and forming R(1) and R(0) ~ where by ¢*(0) = R(1)R(0) 
and P(0) = R(0) +. 


inversion routine. From Lee's work [4 ] the convergence rate appears 


-1 
or by uSing a small sample set of 2N+1 or more 
-] 


Both of these methods require the complicated matrix 


to be satisfactory for any reasonable estimate of 3 *(0) or 4 (0). Ex- 
perience seems to bear this out for the simulation model, The P matrix 
ought to be inversely proportional to time approaching zZéro as 3 
approaches ¢. Phas its most significant effect on the initial rate of 
convergence. Experience indicates that if P is initialized with a 
reasonably large set of values on the main diagonal (~ 10°) , con- 
vergence is satisfactory. Additional work in this area would be useful 
to define an optimal P(0), perhaps in terms of initial convergence rate, 
measurement noise statistics, and identification error as a function of 
Lines 

Convergence tests were run using the scalar observable roll angles 


generated by the fourth order model. A ten second period, damping 


Ne 


coefficient of 0.7 and an excitation constant of twenty were used. The 
sample period was one second. Three cases were investigated using 
the same sample set. 

The first case was with I * suppressed to zeros in all but the last 
element to eliminate numerator dynamics. The convergence of 600 
samples is Shown in fig. 3a,b,c,d for each of the fourth row elements 
of 4 *, The second case was the same as the first but with Gaussian 
additive measurement noise having a standard deviation of one. Ex- 
cellent convergence is shown in fig. 4a,b,c,d. The third case used 
the entire I’ * vector and no meaSurement noise. The residues are 
apparent in fig. 5a,b,c,d. Comparing the eigenvalues of ¢ * with 
those of ¢* with residues indicates a considerable shift, fig. 6. The 
responses to different sample rates in the following section indicates 
that the residues may be minimized sufficiently to provide acceptable 
root locations. 

The true computed values for the plant with a 10 sec period sam- 


pled at 1 sec intervals are as follows, 


0 C 1 0 0 
g*= 0 0 l 0 
0 0 0 i 
[7.1722 JOsGe= =o ety 2, aes 
I * = pl 0s I’ * suppressed = 0 ° 
A245 0 
.900 0 
| 877 | | 877 
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SAMPLE RATE 

The choice of sample rate for discrete stochastic simulation does 
not necessarily follow the supposition that the more rapidly a system 
is sampled the "better". In some cases the existence of an optimal 
sample rate has been shown [ 6 ]. To investigate the accuracy of 
identification by batch processing, a sample set of 500 was generated 
for sample intervals from 0.02 to 4.0 seconds in increments of 0.02 
seconds. The model was excited by a different excitation set for each 
sample rate but with the same statistical properties, normal (0,1) 
"white" noise. The nth row elements of b* for each sample set were 
subtracted from the true ¢* values and plotted versus sample inter- 
val, fig. 7a,b,c,d. For sample intervals less than 0.06 seconds 
R(O) went Singular. Each element appears to have its own best sam- 
ple rate ranging from 1 to 1.5 seconds for a dominant 10 second plant. 
Hence a rule of thumb: Sample at 10 times the dominant frequency! 

To eliminate the effect of varying the excitation set, the above 
procedure was repeated using the same excitation set for each sample 
period, 0.05 to 10 seconds in increments of 0.05 seconds. The re- 
Sulting error curves were the same as illustrated by the 3 *(4,2) curve, 
fig. 8a. The apparent stabilization of the error as the sample period 
approaches the natural period of the system is misleading. The ab- 
solute values of ¢ * become small due to the large damping coefficient 
but percent errors actually increase, 

It is interesting to note that as the dominant period is changed 
the shape of the error curve is the same but shifted in the direction 
of change of the period. The ratios of sample period to natural period 
remain constant at the zero error intercept points. The error curves 
for the 6*(4,2) element for natural periods of 8 and 6 seconds are 
Shown in fig. 8b,c. All of these curves were generated using the full 
I’ * vector and hence contain numerator dynamics and their inherent 
residues. Thus it is pertinent that there is Some optimal sample rate 


that minimizes residue errors, 


SZ 


Errors are not due solely to numerator dynamics. The original 
contention that the choice of sample period is significant may be em- 
phasized by repeating the above runs for the ten second plant with 
suppressed I'*, The resulting identification errors are plotted versus 
Sample period for each element, fig. 9a,b,6e,d. Errors may be sig- 


nificantly different if the excitation, W,, drives the system ata 


ka 
different rate than the sampling (observation) frequency. 
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K-SCALE = 2.0QE+82 UNITS INCH, 
Y-SCALE = 2.086 -@1 UNITS INCH, 


QR af G26 7” G28 810 





Fig. 8. Error = 4 *(4,2) - $ *(4,2) vs Sample Period (0.05 to 10.0 sec) 
using Batch Processing of 500 Samples from the 10, 8, and 6 
sec Dominant Period Plants. One Sequence of 500 Discrete 


Gaussian Values, 0 = 20, was used as Excitation for all Sample 


Periods and all 3 Plants. I * was not Suppressed. (a) Error of 


4 «(4,2) vs Sample Period, 10 sec Plant. (b) Error of $*(4,2) vs 


Sample Period, 8 sec Plant. (c) Error of é *(4,2) vs Sample 
Period, 6 sec Plant. 


(a) Error of 4 * (4,2) vs Sample Period (sec), 10 sec Plant. 
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K-SCALE ~ 2.886-+88 UNITS- INCH. 
4-SCALE = 2.Q8E-81 UNITS INCH 


(b) Error of a *(4,2) vs Sample Period (sec), 8 sec Plant. 
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K-SCALE ~ 2,08€-+88 UNITS~ INCH. 
“-SCALE ~ 2,086-81 UNITS~ INCH. 


(c) Error of 4 *(4,2) vs Sample Period (sec), 6 sec Plant. 
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Error = ¢ * - $ * vs Sample Period (0.05 to 10.0 sec) using 
Batch Processing of 500 Samples from the 10 sec Dominant 
Period Plant with IT * Suppressed (I. * 1,2,3 = 0). Discrete 
Gaussian Excitation, 0 = 20, was applied at the Sample Rate. 
The Same Excitation Set was used for each Different Sample 
Period. (a) Error of r} *(4,1) vs Sample Period. (b) Erropres 

¢@ *(4,2) vs Sample Period. (c) Error of ¢ *(4,3) vs Sample Period. 
(d)’ ETrror or rN *(4,4) vs Sample Period. 


(a) Error of 3 *(4,1) vs Sample Period, “jg sam 


41 


X-SCALE = 2.Q@£+88 UNITS INCH, 
Y-SCALE = 1.88£-@1 UNITS” INCH, 


GO3 


002 


S| 


4D 


O84 806 008 010 





A 
(b) Error of ¢ *(4,2) vs Sample Period, T * 1 202.=_0 


42 


GG1 


001 


003 


49 


G8) QO4 886 008 610 


K-SCALE = 2.00£+88 UNITS INCH. 
‘f-SCALE = 1.086 -91 UNITS/ INCH. 


A 
(c) Error of ¢ *(4,3) vs Sample Period, 1 “| 773) ue 


43 


915 


X-SCALE = 2.00E+80 UNITS~ INCH. 
¥-SCALE = &.Q8E-@2 UNITS~ INCH. 


QAR RR? Q@H G86 088 018 


(d) Error of 3 *(4,4) vs Sample Period, I'*, 2,3 = 0. 
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APPLICATION TO SHIP MOTION ESTIMATION 

As a result of the foregoing work and associated simulation, the 
following proposal is offered as a practical procedure for ship motion 
estimation. 

Initial identification of the order of the pitch and roll plants and 
the corresponding ¢ *'s is to be evaluated for the free dynamic system 
with the ship in calm water. Initial angle may be applied with weights 
or other means. The data thus obtained may be used as a reference and 
to initialize the stochastic routines. The yaw plant may not be iden- 
tified as a free dynamic system since there is no suitable way to apply 
a controlled initial condition which will project on all the eigenvectors 
in the absence of a righting moment. 

Stochastic identification is to be applied whenever the ship is 
underway at sea. Either batch processing or the recursion equations 
may be used. Batch processing should be used at least periodically 
to verify the system order from R(0)~* and to identify Q if required. 

The recursion equations are more Suitable for time sharing and do 
not require storage to accumulate large sample sets of data. Since 
periodic identification requirements are anticipated to follow system 
dynamics, recursive identification may be initiated using the previous 
3% as the best estimate and flagging in P(0). 

The actual estimation of ship attitude follows the identification 
of the ¢ *'s (and Q's). Equation (9) may be applied directly to predict. 
Discrete linear (Kalman) filtering is available. 

The development throughout this paper has assumed that the only 
reliable or convenient measurements available for identification were 
discrete pitch, roll, and yaw angles. If the number of observables 
is increased, the accuracy of the identifications and the estimations 


may be expected to improve [7 ]. 
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CONCLUSIONS 
Identification of ship motion dynamics is feasible by stochastic 
methods, 
The accuracy of the identification is dependent upon the numerator 
dynamics of the true system z-transfer function. The contention 
of Lee [ 4 ] that decorrelation would improve accuracy was not 
true for this work or that of Blackner [7 ]. 
The accuracy is dependent upon the sample period. An optimal 
sample rate exists. 
Identification of an unexcited plant (the free dynamic system) 
having a single initial condition is independent of [T* and hence 
is not affected by numerator dynamics. This should be the source 
of the most reliable identification as a reference, 
Further work should be done using actual ship motion data to 
complete the identification evaluation and to investigate the 
performance of discrete linear (Kalman) filtering as the final 


stage of ship motion estimation. 
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APPENDIX 


CDC FORTRAN 63 Programs and Subroutines for Digital Sim- 


ulation and Identification. 


PHIES?: 


PHIDENT: 


PHICALL: 


PHIDEL: 


PHBE Lass 


PSD: 


REGIP : 


Main program generating the simulation model 
and controlling the identification. 

Recursive identification routine for the scalar 
observable case. 

Batch processing for general observable case, 
Identifies ¢ *, N, and Q. (GAUSS3 is identical 
to RECIP) 

Computes standard form ¢ and I matrices for 
sample data systems. 

Computes ¢ * and I'* matrices for sample 
data systems. 

Computes Power Spectral Density from auto- 
correlation functions, 


Matrix inverSion routine, 


ADD, SUB, TRANS, PROD, PRINTER, and MATREAD: 


SEI 


Convenient subroutines for the indicated 
matrix operations. 
Graphs and prints out 2 dimensional data on 


line printer. 
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